Spatial Ecology and Macroecology

Practical - Week 2

Florencia Grattarola

(Department of Spatial Sciences)

2022-10-03

What are we going to see today?

Gridding biodiversity data, exploring patterns at different resolutions, and making pretty maps.

  1. Gridding biodiversity data
  2. Visualisation

Examples

  • Local: occurrence records 🦌
  • Global: species ranges 🐢

Gridding biodiversity data

Gridding biodiversity data

For all geocomputationgs today we will use the sf package.

install.packages('sf') # install
library(sf) # load
Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE


sf_use_s2(FALSE)
Spherical geometry (s2) switched off

Visualisation

Visualisation

For all the visualisations today we will use the ggplot package. It’s part of the tidyverse.

library(tidyverse)
⁠ggplot(df, aes(x, ⁠

LOCAL Mammal’s of the Czech Republic

Photo of Capreolus capreolus observed in Czech Republic by romanvrbicek licensed CC-BY-NC, via iNaturalist

Mapping occurrence records 🦌

Steps:

  • Get mammal’s of the Czech Republic occurrence records (POINTS)
  • Read the Czech Republic’s grid layer
  • Calculate mammal’s species richness per grid-cell (SR)
  • Visualise SR hotspots in Czech Republic

Mapping occurrence records 🦌

Get the occurrence records

mammalsCZ <- readRDS('data/mammalsCZ.rds')
mammalsCZ
# A tibble: 3,050 × 174
   key    scien…¹ decim…² decim…³ issues datas…⁴ publi…⁵ insta…⁶ publi…⁷ proto…⁸
   <chr>  <chr>     <dbl>   <dbl> <chr>  <chr>   <chr>   <chr>   <chr>   <chr>  
 1 34557… Myodes…    49.8    16.0 cdrou… 50c950… 28eb1a… 997448… US      DWC_AR…
 2 34558… Sus sc…    49.2    16.5 cdrou… 50c950… 28eb1a… 997448… US      DWC_AR…
 3 34558… Myotis…    50.2    16.8 cdrou… 50c950… 28eb1a… 997448… US      DWC_AR…
 4 34558… Lepus …    50.3    15.3 cdrou… 50c950… 28eb1a… 997448… US      DWC_AR…
 5 34558… Capreo…    50.3    15.3 cdrou… 50c950… 28eb1a… 997448… US      DWC_AR…
 6 34561… Cervus…    50.3    14.6 cdrou… 50c950… 28eb1a… 997448… US      DWC_AR…
 7 34562… Capreo…    49.2    17.4 cdrou… 50c950… 28eb1a… 997448… US      DWC_AR…
 8 34564… Dama d…    49.2    16.5 cdrou… 50c950… 28eb1a… 997448… US      DWC_AR…
 9 34564… Lepus …    50.3    15.3 cdrou… 50c950… 28eb1a… 997448… US      DWC_AR…
10 34565… Myocas…    50.0    14.7 cdrou… 50c950… 28eb1a… 997448… US      DWC_AR…
# … with 3,040 more rows, 164 more variables: lastCrawled <chr>,
#   lastParsed <chr>, crawlId <int>, hostingOrganizationKey <chr>,
#   basisOfRecord <chr>, occurrenceStatus <chr>, taxonKey <int>,
#   kingdomKey <int>, phylumKey <int>, classKey <int>, orderKey <int>,
#   familyKey <int>, genusKey <int>, speciesKey <int>, acceptedTaxonKey <int>,
#   acceptedScientificName <chr>, kingdom <chr>, phylum <chr>, order <chr>,
#   family <chr>, genus <chr>, species <chr>, genericName <chr>, …

Let’s keep few fields

mammalsCZ <- mammalsCZ %>% select(species, order, eventDate, 
                                  decimalLongitude, decimalLatitude,
                                  stateProvince, countryCode, datasetName)

Mapping occurrence records 🦌

Let’s transform our data table into an sf object using st_as_sf()

st_as_sf(
  x,
  ...,
  agr = NA_agr_,
  coords,
  wkt,
  dim = "XYZ",
  remove = TRUE,
  na.fail = TRUE,
  sf_column_name = NULL
)

Mapping occurrence records 🦌

Let’s transform our data table into an sf object using st_as_sf()

mammalsCZ %>% 
  filter(!is.na(decimalLongitude) & !is.na(decimalLatitude)) %>% # filter records without coordinates
  st_as_sf(coords=c('decimalLongitude', 'decimalLatitude'),
           crs=4326)
Simple feature collection with 3049 features and 6 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 12.20682 ymin: 48.60141 xmax: 18.79605 ymax: 50.99892
Geodetic CRS:  WGS 84
# A tibble: 3,049 × 7
   species       order event…¹ state…² count…³ datas…⁴            geometry
 * <chr>         <chr> <chr>   <chr>   <chr>   <chr>           <POINT [°]>
 1 Myodes glare… Rode… 2022-0… Pardub… CZ      iNatur… (15.95647 49.80115)
 2 Sus scrofa    Arti… 2022-0… Jihomo… CZ      iNatur… (16.51939 49.19975)
 3 Myotis myotis Chir… 2022-0… Pardub… CZ      iNatur… (16.84201 50.18355)
 4 Lepus europa… Lago… 2022-0… Středo… CZ      iNatur… (15.33439 50.28446)
 5 Capreolus ca… Arti… 2022-0… Středo… CZ      iNatur… (15.33636 50.28574)
 6 Cervus nippon Arti… 2022-0… Středo… CZ      iNatur…  (14.58019 50.3448)
 7 Capreolus ca… Arti… 2022-0… Zlínský CZ      iNatur… (17.36916 49.23076)
 8 Dama dama     Arti… 2022-0… Jihomo… CZ      iNatur…  (16.53508 49.2115)
 9 Lepus europa… Lago… 2022-0… Králov… CZ      iNatur…  (15.2629 50.30738)
10 Myocastor co… Rode… 2022-0… Středo… CZ      iNatur…   (14.653 49.99369)
# … with 3,039 more rows, and abbreviated variable names ¹​eventDate,
#   ²​stateProvince, ³​countryCode, ⁴​datasetName

Note that crs=4326 is related to the EPSG:4326 CRS: WGS 84

Mapping occurrence records 🦌

Let’s transform our data table into an sf object using st_as_sf()

mammalsCZ_sf <- mammalsCZ %>% 
  filter(!is.na(decimalLongitude) & !is.na(decimalLatitude)) %>% # filter records without coordinates
  st_as_sf(coords=c('decimalLongitude', 'decimalLatitude'),
           crs=4326)

Mapping occurrence records 🦌

Get the administrative borders

CZ_borders <- st_read('data/CZE_adm0.shp')
Reading layer `CZE_adm0' from data source 
  `/Users/flograttarola/Documents/GitHub/Spatial-Ecology-and-Macroecology/Practical_classes/Week2_gridding_and plotting/data/CZE_adm0.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 70 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 12.08586 ymin: 48.54084 xmax: 18.86253 ymax: 51.05438
Geodetic CRS:  WGS 84
CZ_borders
Simple feature collection with 1 feature and 70 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 12.08586 ymin: 48.54084 xmax: 18.86253 ymax: 51.05438
Geodetic CRS:  WGS 84
  GADMID ISO     NAME_ENGLI       NAME_ISO       NAME_FAO      NAME_LOCAL
1     59 CZE Czech Republic CZECH REPUBLIC Czech Republic Ceská republika
       NAME_OBSOL NAME_VARIA NAME_NONLA         NAME_FRENC      NAME_SPANI
1 Moravia|Bohemia      Cesko       <NA> République Tchèque República Checa
  NAME_RUSSI         NAME_ARABI NAME_CHINE      WASPARTOF CONTAINS
1      ????? ????????? ????????      ????? Czechoslovakia     <NA>
       SOVEREIGN ISO2  WWW FIPS ISON  VALIDFR VALIDTO AndyID  POP2000     SQKM
1 Czech Republic   CZ <NA>   EZ  203 19930101 Present     65 10271830 78495.16
   POPSQKM      UNREGION1 UNREGION2 DEVELOPING CIS Transition OECD
1 130.8594 Eastern Europe    Europe          2   0          0    1
               WBREGION            WBINCOME        WBDEBT WBOTHER CEEAC CEMAC
1 Europe & Central Asia Upper middle income Less indebted    <NA>     0     0
  CEPLG COMESA EAC ECOWAS IGAD IOC MRU SACU UEMOA UMA PALOP PARTA CACM EurAsEC
1     0      0   0      0    0   0   0    0     0   0     0     0    0       0
  Agadir SAARC ASEAN NAFTA GCC CSN CARICOM EU CAN ACP Landlocked AOSIS SIDS
1      0     0     0     0   0   0       0  1   0   0          1     0    0
  Islands LDC Shape_Leng Shape_Area                       geometry
1       0   0   26.21557   9.813959 MULTIPOLYGON (((13.17648 49...

Mapping occurrence records 🦌

Get the standard grids from the Czech Republic

CZ_grids <- st_read('data/KvadratyCR_JTSK.shp')
Reading layer `KvadratyCR_JTSK' from data source 
  `/Users/flograttarola/Documents/GitHub/Spatial-Ecology-and-Macroecology/Practical_classes/Week2_gridding_and plotting/data/KvadratyCR_JTSK.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 678 features and 10 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -912065.7 ymin: -1229672 xmax: -421278.2 ymax: -928997.8
Projected CRS: S-JTSK / Krovak East North
CZ_grids
Simple feature collection with 678 features and 10 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -912065.7 ymin: -1229672 xmax: -421278.2 ymax: -928997.8
Projected CRS: S-JTSK / Krovak East North
First 10 features:
   OBJECTID ENTITY    AREA PERIMETE KVADRAT         X         Y Shape_Leng
1         1      3 130.032   45.626    4951 -740241.8 -935317.5   45626.41
2         2      5 130.036   45.627    4952 -728665.0 -936923.5   45627.17
3         3      6 130.303   45.675    5051 -741782.9 -946335.6   45675.21
4         4      8 130.040   45.628    4953 -717084.4 -938503.6   45627.96
5         5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
6         6     12 130.311   45.677    5053 -718576.3 -949528.8   45676.59
7         7     14 130.050   45.630    4955 -693912.1 -941586.4   45629.61
8         8     17 130.055   45.630    4956 -682320.5 -943089.1   45630.48
9         9     18 130.319   45.678    5055 -695354.9 -952618.5   45678.09
10       10     19 130.060   45.631    4957 -670725.5 -944566.0   45631.36
   Shape_Area   ID                       geometry
1   130031546 4951 POLYGON ((-745252.4 -928997...
2   130035855 4952 POLYGON ((-733689.7 -930614...
3   130302745 5051 POLYGON ((-746805.7 -940014...
4   130040349 4953 POLYGON ((-722123.3 -932206...
5   130306551 5052 POLYGON ((-735218.5 -941635...
6   130310574 5053 POLYGON ((-723627.4 -943229...
7   130049749 4955 POLYGON ((-698979.2 -935311...
8   130054715 4956 POLYGON ((-687401.8 -936825...
9   130319128 5055 POLYGON ((-700434.2 -946342...
10  130059749 4957 POLYGON ((-675820.7 -938313...

Mapping occurrence records 🦌

Check the layer’s Coordinate Reference System (CRS)

st_crs(CZ_borders)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
st_crs(CZ_grids)
Coordinate Reference System:
  User input: S-JTSK / Krovak East North 
  wkt:
PROJCRS["S-JTSK / Krovak East North",
    BASEGEOGCRS["S-JTSK",
        DATUM["System of the Unified Trigonometrical Cadastral Network",
            ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4156]],
    CONVERSION["Krovak East North (Greenwich)",
        METHOD["Krovak (North Orientated)",
            ID["EPSG",1041]],
        PARAMETER["Latitude of projection centre",49.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8811]],
        PARAMETER["Longitude of origin",24.8333333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8833]],
        PARAMETER["Co-latitude of cone axis",30.2881397527778,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",1036]],
        PARAMETER["Latitude of pseudo standard parallel",78.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8818]],
        PARAMETER["Scale factor on pseudo standard parallel",0.9999,
            SCALEUNIT["unity",1],
            ID["EPSG",8819]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["GIS."],
        AREA["Czechia; Slovakia."],
        BBOX[47.73,12.09,51.06,22.56]],
    ID["EPSG",5514]]

Mapping occurrence records 🦌

We will transform mammalsCZ CRS to S-JTSK / Krovak East North using st_transform()

mammalsCZ_sf <- st_transform(mammalsCZ_sf, crs = st_crs(CZ_grids))
mammalsCZ_sf
Simple feature collection with 3049 features and 6 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -896937 ymin: -1224642 xmax: -435147.4 ymax: -943016.8
Projected CRS: S-JTSK / Krovak East North
# A tibble: 3,049 × 7
   species       order event…¹ state…² count…³ datas…⁴             geometry
 * <chr>         <chr> <chr>   <chr>   <chr>   <chr>            <POINT [m]>
 1 Myodes glare… Rode… 2022-0… Pardub… CZ      iNatur… (-637524.1 -1088504)
 2 Sus scrofa    Arti… 2022-0… Jihomo… CZ      iNatur… (-604525.8 -1159543)
 3 Myotis myotis Chir… 2022-0… Pardub… CZ      iNatur…   (-569727 -1053229)
 4 Lepus europa… Lago… 2022-0… Středo… CZ      iNatur… (-675301.1 -1029784)
 5 Capreolus ca… Arti… 2022-0… Středo… CZ      iNatur… (-675143.9 -1029662)
 6 Cervus nippon Arti… 2022-0… Středo… CZ      iNatur…   (-727699 -1016188)
 7 Capreolus ca… Arti… 2022-0… Zlínský CZ      iNatur…   (-542595 -1162493)
 8 Dama dama     Arti… 2022-0… Jihomo… CZ      iNatur… (-603247.5 -1158367)
 9 Lepus europa… Lago… 2022-0… Králov… CZ      iNatur… (-680037.4 -1026620)
10 Myocastor co… Rode… 2022-0… Středo… CZ      iNatur… (-727767.4 -1055585)
# … with 3,039 more rows, and abbreviated variable names ¹​eventDate,
#   ²​stateProvince, ³​countryCode, ⁴​datasetName

Mapping occurrence records 🦌

Let’s plot the sf objects

ggplot() + 
  geom_sf(data=CZ_borders, fill='white') +
  geom_sf(data=CZ_grids, fill=NA) + 
  geom_sf(data=mammalsCZ_sf)

Mapping occurrence records 🦌

To calculate number of records (N) and number of species per grid-cell (SR), we will use st_join()

st_join(
  x,    #   object of class sf
  y,    #   object of class sf
  join = st_intersects,
  ...,
  suffix = c(".x", ".y"),
  left = TRUE,
  largest = FALSE
)

Mapping occurrence records 🦌

To calculate number of records (N) and number of species per grid-cell (SR), we will use st_join().

st_join(CZ_grids,       # GRID
        mammalsCZ_sf)   # POINTS
Simple feature collection with 3350 features and 16 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -912065.7 ymin: -1229672 xmax: -421278.2 ymax: -928997.8
Projected CRS: S-JTSK / Krovak East North
First 10 features:
    OBJECTID ENTITY    AREA PERIMETE KVADRAT         X         Y Shape_Leng
1          1      3 130.032   45.626    4951 -740241.8 -935317.5   45626.41
2          2      5 130.036   45.627    4952 -728665.0 -936923.5   45627.17
3          3      6 130.303   45.675    5051 -741782.9 -946335.6   45675.21
4          4      8 130.040   45.628    4953 -717084.4 -938503.6   45627.96
5          5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
5.1        5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
5.2        5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
5.3        5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
5.4        5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
5.5        5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
    Shape_Area   ID             species        order           eventDate
1    130031546 4951                <NA>         <NA>                <NA>
2    130035855 4952                <NA>         <NA>                <NA>
3    130302745 5051                <NA>         <NA>                <NA>
4    130040349 4953                <NA>         <NA>                <NA>
5    130306551 5052 Capreolus capreolus Artiodactyla 2020-06-23T00:00:00
5.1  130306551 5052 Capreolus capreolus Artiodactyla 2020-06-29T00:00:00
5.2  130306551 5052 Capreolus capreolus Artiodactyla 2020-06-27T00:00:00
5.3  130306551 5052     Lepus europaeus   Lagomorpha 2020-06-29T00:00:00
5.4  130306551 5052 Capreolus capreolus Artiodactyla 2020-06-24T00:00:00
5.5  130306551 5052     Lepus europaeus   Lagomorpha 2020-06-26T00:00:00
    stateProvince countryCode datasetName                       geometry
1            <NA>        <NA>        <NA> POLYGON ((-745252.4 -928997...
2            <NA>        <NA>        <NA> POLYGON ((-733689.7 -930614...
3            <NA>        <NA>        <NA> POLYGON ((-746805.7 -940014...
4            <NA>        <NA>        <NA> POLYGON ((-722123.3 -932206...
5            <NA>          CZ        <NA> POLYGON ((-735218.5 -941635...
5.1          <NA>          CZ        <NA> POLYGON ((-735218.5 -941635...
5.2          <NA>          CZ        <NA> POLYGON ((-735218.5 -941635...
5.3          <NA>          CZ        <NA> POLYGON ((-735218.5 -941635...
5.4          <NA>          CZ        <NA> POLYGON ((-735218.5 -941635...
5.5          <NA>          CZ        <NA> POLYGON ((-735218.5 -941635...

Mapping occurrence records 🦌

Now, we need to summarise the data per grid-cell.
Luckily, Tidyverse methods also work for sf objects :)


We will do it with group_by() and summarise().

First, we group by grid-cells and then we count values per grid.

Mapping occurrence records 🦌

Let’s calculate number of records (N) and number of species per grid-cell (SR)

st_join(CZ_grids, mammalsCZ_sf) %>% 
  group_by(OBJECTID) %>% 
  summarise(N=n(),
            SR=n_distinct(species))
Simple feature collection with 678 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -912065.7 ymin: -1229672 xmax: -421278.2 ymax: -928997.8
Projected CRS: S-JTSK / Krovak East North
# A tibble: 678 × 4
   OBJECTID     N    SR                                                 geometry
      <int> <int> <int>                                            <POLYGON [m]>
 1        1     1     1 ((-745252.4 -928997.8, -733689.7 -930614.9, -733689.7 -…
 2        2     1     1 ((-733689.7 -930614.9, -722123.3 -932206.2, -722123.3 -…
 3        3     1     1 ((-746805.7 -940014.3, -735218.5 -941635, -735218.5 -94…
 4        4     1     1 ((-722123.3 -932206.2, -710553.1 -933771.6, -710553.1 -…
 5        5    27     3 ((-735218.5 -941635, -723627.4 -943229.9, -723627.4 -94…
 6        6     4     3 ((-723627.4 -943229.9, -712032.7 -944798.9, -712032.7 -…
 7        7     1     1 ((-698979.2 -935311.3, -687401.8 -936825.2, -687401.8 -…
 8        8     1     1 ((-687401.8 -936825.2, -675820.7 -938313.3, -675820.7 -…
 9        9     1     1 ((-700434.2 -946342, -688832.2 -947859.3, -688832.2 -94…
10       10     1     1 ((-675820.7 -938313.3, -664236.1 -939775.7, -664236.1 -…
# … with 668 more rows

Mapping occurrence records 🦌

Let’s store the output into a new object CZ_mammals_grids.

mammalsCZ_grids <- st_join(CZ_grids, mammalsCZ_sf) %>%
  group_by(OBJECTID) %>%
  summarise(N=ifelse(n_distinct(species, na.rm = TRUE)==0, 0, n()),
            SR=n_distinct(species, na.rm = TRUE))

Mapping occurrence records 🦌

Finally, let’s plot this.
We will do it using geom_sf() a ggplot function to visualise sf objects.

Mapping occurrence records 🦌

ggplot() + 
  geom_sf(data=mammalsCZ_grids) +
  coord_sf(crs=4326)

Where are the nice colours?

Mapping occurrence records 🦌

We need to indicate which column from the object should the grids be filled with.

Mapping occurrence records 🦌

ggplot() + 
  geom_sf(data=mammalsCZ_grids, aes(fill=SR)) +
  coord_sf(crs=4326)

Mapping occurrence records 🦌

Let’s get a nicer color scale.

Bare in mind that we see colors differently. Thus, it’s important to consider colorblind safe palettes

Check out The R Graph Gallery for many cool graphic options with ggplot

Mapping occurrence records 🦌

ggplot() + 
  geom_sf(data=mammalsCZ_grids, aes(fill=SR)) +
  scale_fill_fermenter(palette ='YlOrBr', n.breaks=9, direction = 1)+
  geom_sf(data=CZ_borders, fill=NA) +
  coord_sf(crs=4326) +
  theme_bw()

Now, here’s a better plot :)

Mapping occurrence records 🦌

The hotspots are located in cities, how can they be the richest areas for mammals?

Let’s have a look at the sampling effort (N: number of records) and compare both layers.

Mapping occurrence records 🦌

ggplot() + 
  geom_sf(data=mammalsCZ_grids, aes(fill=N)) +
  scale_fill_fermenter(palette ='YlGnBu', n.breaks=9, direction = 1) +
  geom_sf(data=CZ_borders, fill=NA) +
  coord_sf(crs=4326) +
  theme_bw()

Mapping occurrence records 🦌

What can you say about the hotspots of species richness we found?

Mapping occurrence records 🦌

  • Get mammal’s of the Czech Republic occurrence records (POINTS)
  • Get the Czech Republic’s grid
  • Calculate mammal’s species richness per grid-cell (SR)
  • Visualise SR hotspots in Czech Republic


We are done! No it’s your turn :)

GLOBAL Testudines of the World

Photo of Emys orbicularis observed in France by Martin Costechareire licensed CC-BY-NC, via iNaturalist

Mapping species ranges 🐢

Steps:

  • Read testudines’ IUCN range maps (POLYGONS)
  • Get a world map
  • Create a global 1-degree grid
  • Calculate testudines’s species richness per grid-cell (SR)
  • Visualise SR global hotspots

Mapping species ranges 🐢

Read testudines’ IUCN range maps (POLYGONS)

testudines <- st_read('data/testudines/data_0.shp')
Reading layer `data_0' from data source 
  `/Users/flograttarola/Documents/GitHub/Spatial-Ecology-and-Macroecology/Practical_classes/Week2_gridding_and plotting/data/testudines/data_0.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 169 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -179.999 ymin: -47.55406 xmax: 179.999 ymax: 61.50515
Geodetic CRS:  WGS 84


Polygons were downloaded from IUCN Spatial Data Download’s page.

Mapping species ranges 🐢

What do the data look like?

testudines
Simple feature collection with 169 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -179.999 ymin: -47.55406 xmax: 179.999 ymax: 61.50515
Geodetic CRS:  WGS 84
First 10 features:
   ASSESSMENT ID_NO                 BINOMIAL PRESENCE ORIGIN SEASONAL
1      495907 10041      Heosemys annandalii        1      1        1
2      499158 10825    Indotestudo forstenii        1      1        1
3      499618 10950    Pangshura sylhetensis        1      1        1
4      505402 12124         Lissemys scutata        1      1        1
5      507698 12695      Malaclemys terrapin        1      1        1
6      507698 12695      Malaclemys terrapin        1      5        1
7      508210 12696   Malacochersus tornieri        1      1        1
8      508518 12775        Manouria impressa        1      1        1
9      511526 13038 Melanochelys tricarinata        1      1        1
10     511526 13038 Melanochelys tricarinata        6      1        1
        COMPILER YEAR                      CITATION
1  Anders Rhodin 2020 Chelonian Research Foundation
2  Anders Rhodin 2020 Chelonian Research Foundation
3  Anders Rhodin 2020 Chelonian Research Foundation
4  Anders Rhodin 2020 Chelonian Research Foundation
5  Anders Rhodin 2018 Chelonian Research Foundation
6  Anders Rhodin 2018 Chelonian Research Foundation
7  Anders Rhodin 2018 Chelonian Research Foundation
8  Anders Rhodin 2020 Chelonian Research Foundation
9  Anders Rhodin 2018 Chelonian Research Foundation
10 Anders Rhodin 2018 Chelonian Research Foundation
                                 LEGEND SUBSPECIES SUBPOP DIST_COMM ISLAND
1                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
2                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
3                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
4                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
5                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
6  Extant & Origin Uncertain (resident)       <NA>   <NA>      <NA>   <NA>
7                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
8                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
9                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
10                   Presence Uncertain       <NA>   <NA>      <NA>   <NA>
   TAX_COMM                       geometry
1      <NA> MULTIPOLYGON (((102.6212 5....
2      <NA> MULTIPOLYGON (((120.8806 1....
3      <NA> MULTIPOLYGON (((92.29536 22...
4      <NA> MULTIPOLYGON (((95.29987 15...
5      <NA> MULTIPOLYGON (((-81.45417 2...
6      <NA> MULTIPOLYGON (((-64.7296 32...
7      <NA> MULTIPOLYGON (((32.93152 -9...
8      <NA> MULTIPOLYGON (((101.684 5.0...
9      <NA> MULTIPOLYGON (((85.58097 25...
10     <NA> MULTIPOLYGON (((92 21.59583...

Mapping species ranges 🐢

Now, let’s get a world map. For this we will use the package rnaturalearth.

install.packages(c('rnaturalearth', 'rnaturalearthdata'))
library(rnaturalearth)
library(rnaturalearthdata)

Mapping species ranges 🐢

We will download a shapefile of the world at medium scale and we will combine all countries into one unique polygon.

world <- ne_countries(scale = 50, returnclass='sf')
world <- st_union(world) %>% st_make_valid() %>% st_cast()


Let’s see how it looks

Mapping species ranges 🐢

ggplot() + 
  geom_sf(data=world, fill='white')

Mapping species ranges 🐢

Before working with this map, we need to project the layer: from lat/lon to equal area projection.


Earth is not flat :) Projections help us represent the two-dimensional curved surface of the globe into 2D space. There are many ways to do this (map projections).

Mapping species ranges 🐢

Equal-area maps preserve area measure, generally distorting shapes in order to do that.

Mapping species ranges 🐢

We will choose Equal Earth (EPSG:8857).

Mapping species ranges 🐢

Let’s transform the data (both the world and the testudines’ layers)

world_ea <- st_transform(world, crs = 'EPSG:8857') %>% 
  st_make_valid %>% st_cast
testudines_ea <- st_transform(testudines, crs = 'EPSG:8857') %>% 
  st_make_valid %>% st_cast


And double check everything’s alright

st_crs(world_ea, parameters = TRUE)$epsg
[1] 8857
st_crs(world_ea, parameters = TRUE)$ud_unit
1 [m]
st_crs(testudines_ea, parameters = TRUE)$epsg
[1] 8857
st_crs(testudines_ea, parameters = TRUE)$ud_unit
1 [m]

Mapping species ranges 🐢

Now, let’s create 1-degree grid-cells for the entire world.
We will do this using the function st_make_grid()

st_make_grid(
  x,
  cellsize = c(diff(st_bbox(x)[c(1, 3)]), diff(st_bbox(x)[c(2, 4)]))/n,
  offset = st_bbox(x)[c("xmin", "ymin")],
  n = c(10, 10),
  crs = if (missing(x)) NA_crs_ else st_crs(x),
  what = "polygons",
  square = TRUE,
  flat_topped = FALSE
)

Mapping species ranges 🐢

For the cellsize we will chose 1 degree, which is ~100km (=100,000m)

world_grid <- st_make_grid(world_ea, 
                           cellsize=100000,
                           square=FALSE) %>%  # this will make hexagons
  st_intersection(world_ea) %>% 
  st_sf(gridID=1:length(.))


We have the grid, we are ready to calculate richness.

Mapping species ranges 🐢

Mapping species ranges 🐢

To make things faster, let’s get a smaller dataset of Testudines

testu <- testudines_ea %>% 
  filter(LEGEND=='Extant (resident)') %>% 
  group_by(ID_NO) %>% 
  summarise(BINOMIAL=unique(BINOMIAL),
            PRESENCE=sum(PRESENCE)) %>% 
  ungroup() %>% st_cast()


This can take a while

Mapping species ranges 🐢

Let’s plot 10 species as an example

Mapping species ranges 🐢

To calculate number of records (N) and number of species per grid-cell (SR), we will use st_join().

testudines_grid <- st_join(world_grid, testu) %>%
  group_by(gridID) %>%
  summarise(N=ifelse(n_distinct(BINOMIAL, na.rm = TRUE)==0, 0, n()),
            SR=n_distinct(BINOMIAL, na.rm = TRUE))


This will take a while


And now, we plot!

Mapping species ranges 🐢

Here’s the results :)

testudines_grid
Simple feature collection with 20146 features and 3 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: -16920570 ymin: -8392928 xmax: 16920570 ymax: 8315130
Projected CRS: WGS 84 / Equal Earth Greenwich
# A tibble: 20,146 × 4
   gridID     N    SR                                                          .
    <int> <dbl> <int>                                             <GEOMETRY [m]>
 1      1     0     0                                 POINT (-16920565 -2061608)
 2      2     2     2 MULTIPOLYGON (((-16903934 -2109412, -16904928 -2110323, -…
 3      3     2     2 MULTIPOLYGON (((-16905520 -2108496, -16907015 -2103886, -…
 4      4     2     2 POLYGON ((-16782581 -2193350, -16785819 -2191050, -167869…
 5      5     3     3 POLYGON ((-16805308 -1830105, -16810425 -1830932, -168160…
 6      6     2     2 POLYGON ((-16783058 -2419039, -16784544 -2418515, -167857…
 7      7     2     2 MULTIPOLYGON (((-16779389 -2191507, -16775928 -2194312, -…
 8      8     2     2 MULTIPOLYGON (((-16749600 -2289778, -16749039 -2290664, -…
 9      9     2     2 POLYGON ((-16692191 -601861.4, -16691700 -601711.1, -1669…
10     10     2     2 POLYGON ((-16658328 -2432319, -16657130 -2432783, -166574…
# … with 20,136 more rows

Mapping species ranges 🐢

References

Any doubts?